%% Figure2.m
%  This script reproduces figure 2.B and some statistical analyses,
%  computes the correlation coeffcients and displays computational times
%
% Estimated run time: 10 seconds.
clc
close all
clear
load DataSet.mat

for chromo=1:size(Unici_sani,1)
clear CR CRS CR_event CR_ordered CR_ordered2 ISI Rep_Duration i count nrep
cd(['.\SAVED_FILES_Generation_1_Chromosome_' num2str(chromo-1) ])
CR=load('CR.dat');
count=1;
nrep=70;
Rep_Duration=0.8;
repduration=800;
ISI=600;
if numel(CR)>0   
CR_ordered(1,1)=CR(1);
for i=2:size(CR,1)
    if CR(i)~=CR(i-1)
        CR_ordered(count+1,1)=CR(i);
        count=count+1;
    end
end

% MOBILE
count=1;
for i=1:nrep
    CR_ordered2=CR_ordered(find(CR_ordered<=(i*Rep_Duration)));
    CRS=CR_ordered2(find(CR_ordered2>(i*Rep_Duration-Rep_Duration*10)));
    CR_event(count)=size(CRS,1)*10;
    count=count+1;
   
end
    CR_sani(:,chromo)=CR_event;
else
    CR_sani(:,chromo)=zeros(nrep,1);
end

clear i j k total filename count conta
cd('..')
end
clear CR CRS CR_event CR_ordered CR_ordered2 ISI Rep_Duration i count nrep chromo

CR_sani_num=[];
for chromo=1:size(Unici_sani,1)
    volte=Numerosita_sani(chromo);
    CR_sani_num=[CR_sani_num;repmat(CR_sani(:,chromo)',volte,1)];
end
clear volte chromo

figureFullScreen
stdshade_median(CR_sani_num,0.2,'b')
hold on
stdshade_median([MobileM1 MobileM3]',0.2,'k')
title('figure 2.B - session_1','fontsize',20)
set(gca,'fontsize',20);
box off
axis([1 70 0 100 ])
set(gca,'fontsize',20)
legend('Model Data','','Experimental Data','Location','Best')
xlabel('Trial','fontsize',20)

%% SANI2 

for chromo=1:size(Unici_sani2,1)
clear CR CRS CR_event CR_ordered CR_ordered2 ISI Rep_Duration i count nrep
cd(['.\SAVED_FILES_Generation_2_Chromosome_' num2str(chromo-1) ])
CR=load('CR.dat');
count=1;
nrep=70;
Rep_Duration=0.8;
ISI=600;
if numel(CR)>0   
CR_ordered(1,1)=CR(1);
for i=2:size(CR,1)
    if CR(i)~=CR(i-1)
        CR_ordered(count+1,1)=CR(i);
        count=count+1;
    end
end

% MOBILE
count=1;
for i=1:nrep
    CR_ordered2=CR_ordered(find(CR_ordered<=(i*Rep_Duration)));
    CRS=CR_ordered2(find(CR_ordered2>(i*Rep_Duration-Rep_Duration*10)));
    CR_event(count)=size(CRS,1)*10;
    count=count+1;
   
end
    CR_sani2(:,chromo)=CR_event;
else
    CR_sani2(:,chromo)=zeros(nrep,1);
end

clear i j k total filename count conta
cd '..'
end
clear CR CRS CR_event CR_ordered CR_ordered2 ISI Rep_Duration i count nrep chromo

CR_sani_num2=[];
for chromo=1:size(Unici_sani2,1)
    volte=Numerosita_sani2(chromo);
    CR_sani_num2=[CR_sani_num2;repmat(CR_sani2(:,chromo)',volte,1)];
end
clear volte chromo
figureFullScreen
stdshade_median(CR_sani_num2,0.2,'g')
hold on
stdshade_median(MobileM2',0.2,'k')
title('figure 2.B - session_{2 sham}','fontsize',20)
set(gca,'fontsize',20);
box off
axis([1 70 0 100 ])
set(gca,'fontsize',20)
legend('Model Data','','Experimental Data','Location','Best')
xlabel('Trial','fontsize',20)

%% TMS2

for chromo=1:size(Unici_tms2,1)
clear CR CRS CR_event CR_ordered CR_ordered2 ISI Rep_Duration i count nrep
cd(['.\SAVED_FILES_Generation_3_Chromosome_' num2str(chromo-1) ])
CR=load('CR.dat');
count=1;
nrep=70;
Rep_Duration=0.8;
ISI=600;
if numel(CR)>0   
CR_ordered(1,1)=CR(1);
for i=2:size(CR,1)
    if CR(i)~=CR(i-1)
        CR_ordered(count+1,1)=CR(i);
        count=count+1;
    end
end

% MOBILE
count=1;
for i=1:nrep
    CR_ordered2=CR_ordered(find(CR_ordered<=(i*Rep_Duration)));
    CRS=CR_ordered2(find(CR_ordered2>(i*Rep_Duration-Rep_Duration*10)));
    CR_event(count)=size(CRS,1)*10;
    count=count+1;
   
end
    CR_tms2(:,chromo)=CR_event;
else
    CR_tms2(:,chromo)=zeros(nrep,1);
end

clear i j k total filename count conta
cd '..'
end
clear CR CRS CR_event CR_ordered CR_ordered2 ISI Rep_Duration i count nrep chromo

CR_tms_num2=[];
for chromo=1:size(Unici_tms2,1)
    volte=Numerosita_tms2(chromo);
    CR_tms_num2=[CR_tms_num2;repmat(CR_tms2(:,chromo)',volte,1)];
end
clear volte chromo
figureFullScreen
stdshade_median(CR_tms_num2,0.2,'r')
hold on
stdshade_median(MobileM4',0.2,'k')
title('figure 2.B - session_{2 tbs}','fontsize',20)
set(gca,'fontsize',20);
box off
axis([1 70 0 100 ])
set(gca,'fontsize',20)
legend('Model Data','','Experimental Data','Location','Best')
xlabel('Trial','fontsize',20)

%% DISPLAY PEARSON CORRELATION COEFFIFIENTS
%  Between experimental and model medians

rho1=corr2(median(CR_sani_num)',median([MobileM1 MobileM3]')');
rho2=corr2(median(CR_sani_num2)',median(MobileM2')');
rho3=corr2(median(CR_tms_num2)',median(MobileM4')');

disp('%%% Pearson Correlation Indexes %%%')
disp(['Correlation Session1 = ' num2str(rho1)])
disp(['Correlation Session2 sham = ' num2str(rho2)])
disp(['Correlation Session2 tbs = ' num2str(rho3)])


%% DISPLAY COMPUTATIONAL TIME

gen_num = 1;
for chromo=1:size(Unici_sani,1)
    cd(['.\SAVED_FILES_Generation_' num2str(gen_num) '_Chromosome_' num2str(chromo-1)])
    filename = 'log.dat';
    delimiter = '\t';
    startRow = 13;
    endRow = 13;
    formatSpec = '%s%*s%*s%*s%*s%*s%[^\n\r]';
    fileID = fopen(filename,'r');
    dataArray = textscan(fileID, formatSpec, endRow-startRow+1, 'Delimiter', delimiter, 'HeaderLines', startRow-1, 'ReturnOnError', false);
    fclose(fileID);
    clearvars filename delimiter startRow endRow formatSpec fileID;
    tempo=char(dataArray{1});
    i=0;
    while isempty(str2num(tempo(15:20-i)))
        i=i+1;
    end
    tempo2(gen_num,chromo)=str2num(tempo(15:20-i));
    cd('..')
end

gen_num = 2;
for chromo=1:size(Unici_sani2,1)
    cd(['.\SAVED_FILES_Generation_' num2str(gen_num) '_Chromosome_' num2str(chromo-1)])
    filename = 'log.dat';
    delimiter = '\t';
    startRow = 13;
    endRow = 13;
    formatSpec = '%s%*s%*s%*s%*s%*s%[^\n\r]';
    fileID = fopen(filename,'r');
    dataArray = textscan(fileID, formatSpec, endRow-startRow+1, 'Delimiter', delimiter, 'HeaderLines', startRow-1, 'ReturnOnError', false);
    fclose(fileID);
    clearvars filename delimiter startRow endRow formatSpec fileID;
    tempo=char(dataArray{1});
    i=0;
    while isempty(str2num(tempo(15:20-i)))
        i=i+1;
    end
    tempo2(gen_num,chromo)=str2num(tempo(15:20-i));
    cd('..')
end

gen_num = 3;
for chromo=1:size(Unici_tms2,1)
    cd(['.\SAVED_FILES_Generation_' num2str(gen_num) '_Chromosome_' num2str(chromo-1)])
    filename = 'log.dat';
    delimiter = '\t';
    startRow = 13;
    endRow = 13;
    formatSpec = '%s%*s%*s%*s%*s%*s%[^\n\r]';
    fileID = fopen(filename,'r');
    dataArray = textscan(fileID, formatSpec, endRow-startRow+1, 'Delimiter', delimiter, 'HeaderLines', startRow-1, 'ReturnOnError', false);
    fclose(fileID);
    clearvars filename delimiter startRow endRow formatSpec fileID;
    tempo=char(dataArray{1});
    i=0;
    while isempty(str2num(tempo(15:20-i)))
        i=i+1;
    end
    tempo2(gen_num,chromo)=str2num(tempo(15:20-i));
    cd('..')
end

tempo2(find(tempo2==0))=NaN;

disp('%%% Computational Times %%%')
disp(['Session1 computational time: ' num2str(nanmean(tempo2(1,:))) ' � ' num2str(nanstd(tempo2(1,:))) 'seconds'])
disp(['Session2 sham computational time: ' num2str(nanmean(tempo2(2,:))) ' � ' num2str(nanstd(tempo2(2,:))) 'seconds'])
disp(['Session2 tbs computational time: ' num2str(nanmean(tempo2(3,:))) ' � ' num2str(nanstd(tempo2(3,:))) 'seconds'])
